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ABSTRACT 

Galaxies can form in a sufficiently deep gravitational potential so that efficient 
gas cooling occurs. We estimate that such potential is provided by a halo of mass 
M > Mc Si 7.0 X 10'^ (Ar(z)(l +zy'T^''^ h~^ Mq, where Adz) is the mean overdensity of 
spherically virialized objects formed at redshift z, and Mc ~ 4.0 x 10^' h~^ Mq at z = 0. 
Based on this criterion, our galaxy samples are constructed from cosmology simulation 
data by using HiFOF to select subhalos in those FOF halos that are more massive than 
Mc. There are far more dark subhalos than galaxy-hosting subhalos. Several tests 
against observations have been performed to examine our galaxy samples, including 
the diff"erential galaxy mass functions, the galaxy space density, the projected two 
point correlation functions (CF), the HODs, and the kinematic pair fractions. These 
tests show good agreements. Based on the consistency with observations, our galaxy 
sample is believed to correctly represent galaxies in real universe, and can be used to 
study other unexplored galaxy properties. 

Subject headings: cosmology: theory — galaxies: formation — galaxies: halos — 
large-scale structure of universe — galaxies: evolution — methods: numerical 
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1. Introduction 



The standard theoretical model for galaxy formation basically consists of two main elements, 
the cold dark matter (CDM) model and a dark energy field (which may take the form of a cos- 
mological constant, A). The fundamental assumption made in the theory is that structure grew 
from weak density fluctuations present in the otherwise homogeneous and rapidly expanding early 
Universe. Afterwards, the gravitational instability drives these fluctuations into nonlinear regime 
in a bottom up fashion and they gradually grow to a wealth of structures to day. Recent divers e 
cosmological studies, no m atter from large-s cale structure observations (e.g., ISpergel et al.ll2003b 



from supernova data (e.g., iRiess et al.lll998h . or from light element abundance (e.g., IWhite et al 



1993b . all seem to support the concordance model. Due to the highly nonlinear nature in the col- 
lapse of fluctuations and the subsequent hierarchical build-up of structure, numerical simulations 
come to play an indispensable role. 

In the past thirty years, galaxy clustering has been intensively studied, as the clustering of 
galaxies has long been an essential testing ground for various cosmological models and galaxy for- 
mation scenarios. Especially due to the advent of modem computers, high resolution simulations 
become possible. Cosmological N-body simulations have developed into a powerful tool for calcu- 
lating the gravitational clustering of coUisionless dark matter from specified initial conditions. As 
the resolution of sim ulations increases, some numerical problems such as the overmerging problem 
(|Klypin et al.lll999|) can be overcome, and the galactic size scale in a large scale structure simula- 
tion can be resolvable. In addition, the growth of galaxy redshift surveys has led to measurements 
of increasing precision and detail. All these works have made the comparisons between theories 
and observations possible. 

There are four main approaches used in cosmological simulations to study the galaxy clus- 
tering and properties. The fi rst one is N-body plu s hydrodynamical simulation including gas and 
dark matter particles (e.g., [Weinberg et al.l 120041) . The second method is a hybrid method that 
combines N-body simulations o f the dark matter com ponent with semi-analytic treatments of the 
galaxy formation physics (e.g., Springel et al.ll2005|) . Another approach is using high-resolution 
coUisionless N-body simulations that identify galaxies with "subhalos" in the dark matter distri- 
bution (e.g., iKly pin et al. 1999 : Kravtsov et al. 2004; Conroy et al. 2006). The last one is the so 
called Halo Occupation Distribution (HOD) approach wh ich gives a purely statistical description of 
how dark matter h alos are populated with galaxies (e.g., iBerlind et al.ll2003l: iKravtsov et al.ll2004l : 
Zheng et aPboOS ). It is known that the cold dark matter is the dominant mass component and it in- 
teracts only through gravity. To some extent, gravita tional dynamics alone should explain the basic 
features of galaxy clustering. iNeyrinck et al.l (|2004|) tried to understand how the infrared-selected 
galaxies populate dark matter halos, paying special attention to the method of halo identification 
in simulations. They tested the hypothesis that baryonic physics negligibly aff"ects the distribu- 
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tion of subhalos down to the smallest scales yet obs erved and suc c essful ly reproduced the Point 
Source Catalogue Redshift (PSCz) power spectrum. iBerlind et alj (|2003[) reported that the HOD 
prediction of two methods, a semi-analytic model and gasdynamics simulations, agree remarkably 
well for samples of the same space density. This result indirectly supported the idea that the HOD, 
and hence galaxy clustering, is driven primarily by gra vitational dynamics ra ther than by processes 



such as cooling and star formation. On the other hand. lKravtsov et al 



bound density maxima (BDM) halo-finding algorithm (IKlypin et al 



(120041) adopt a variant of the 



19991) to identify halos and 



subhalos and use the maximum circular velocity V,nax as a proxy of halo mass for selecting galaxy 
sample. Instead of selecting objects in a given range of V,nax, at each epoch they selected objects 
of a set of number densities consistent with observational number densities and corresponding to 
(redshift dependent) thresholds in maximum circular velocity, i.e., n,(> V,„ax)- Their results show 
that the dependence of correlation amplitude on the galaxy number density in their sarn ple is in 
general agreement with results from the Sloan Digital Sky Survey. IConroy et al.l (|2006h instead 
use the maximum circular velocity at the time of accretion, V^^^^, for subhalos, and the results 
show agreement with the observed galaxy clustering in the SDSS data at z « and in the DEEP2 
samples at z « 1 over the range of separations, 0.1 < rp/{h~^Mpc) < 10.0. However, this work 
lacks the information of the subhalo mass. Motivated by the results of these previous works, an 
approach similar to the galaxy identification with subhalos is adopted in this study. Despite some 
of our simulations include gas particles, we shall ignore the gas component in the identification of 
galaxies. Our approach in fact combines a subhalo finding algorithm, HiFOF, and a galaxy forma- 
tion model. It 'directly' locates where the galaxy should be formed and shows agreements with 
various observation results, such as galaxy mass function, two-point correlation function, etc. This 
result indicates that our galaxy identification model not only can successfully locate subhalos in 
simulations that reside inside host halos in real universe but also can give a galaxy mass function 
close to the observational one. 

The paper is structured as follows. In section 2, we describe the simulation and the halo- 
finding algorithm in use. Thereafter, the galaxy halo model we adopt is discussed. In Section 
3, we compare our galaxy samples with observations including the differential mass functions, 
the number density evolution, the correlation functions at different redshift, the halo occupation 
distribution, and the pair fraction. Finally, we discuss and summarize our results in Section 4. 
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2. Theoretical models 



2.1. Simulations 



Our simulations have been evolved in the concordance flat ACDM model: = 0-3, = 0.7 
Q.h = 0.05, and h = 0.66, where Qq, Q.^, Qi, are the present-day matter, vacuum, and baryon den- 
sities and h is the dimensionless Hubble constant defined as Hq = 100 h km Mpc~^. These pa 



rameters are consistent with current observations on cosmol ogical parameters (|Spergel et al.ll2003 
Tegmark et all2004l : lseljak et alboosl : Isanchez et al.ll2006l) with Q^^ = 0.20 + 0.020, Q,, = 0.042 
+ 0.002, Q.A = 0.76 + 0.020, h = 0.74 + 0.02. Four realizations with different box sizes and particle 
numbers were analyzed and compare d to investigate the bo undary and resolu tion effects as w ell as 
to have better statistics. GADGET 1 dSpringel et allboOlbl) and GADGET2 (lspringel||2005|) were 
both employed to conduct the simulations. Our simulations all started at redshift z = 100 and 
evolved to z = 0, and their density spectra were normalized to erg = 0.94, where cr^ is the rms 
fluctuation in sphere of 8 Mpc comoving radius. The first simulation (ACDMiooo) follows the 
evolution of 256^ dark matter particles and 256^ gas particles in a 100 Mpc box on a side. The 
mass of a dark matter particle is ntdm = 4.125 xlO^ Mq, while the mass of a gas particle is 
nigas = 8.25 X 10^ Mq. In addition. We adopted a softening length switching from a comoving 
scale to a physical scale at z = 2.3. The softening length e was 20 kpc (comoving) before red- 
shift z = 2.3. After that, e was switched to 6 kpc (physical). Thus, the highest force resolution 
is 6 h'^ kpc. The second simulation (ACDM200) follows the evolution of 512^ dark matter particles 
and 512^ gas particles in the same cosmology but in a 200 Mpc box on a side. Both dark matter 
particles and gas particles have the same masses respectively as in the first simulation, but in this 
run e keeps constant at 10 h'^ kpc (comoving). The third simulation (ACDMioo^) evolves 512^ 
pure dark matter particles in a 100 Mpc box on a side. The mass of a dark matter particle is 
nidm = 6.188 xlO*^ h'^ Mq. We also adopt a softening length switching scheme in this simulation. 
However, the softening length e was set to be 10 kpc (comoving) before redshift z = 2.3 and 
was modified to 3 h'^ kpc (physical) thereafter. The final simulation (ACDMiooc) was run with 
512^ dark matter particles and 512^ gas particles in a 100 Mpc box on a side, and m^,,, = 5.156 
xlO^ h^'^ Mq, nigas = 1.031 xlO^ h'^ Mq, and e was kept constant at 10 h'^ kpc (comoving). The 
parameters of our simulations are summarized in Table [IJ 



2.2. Hierarchical Friends-of-Friends Algorithm (HiFOF) 



There are many widely used algorith ms for id entifying the subst ructures, such as Bound Den- 
sity M axima (BDM ; iKlvpi n et al.lll999h . SKID dstadel et all 1 19971) . and subfind d Springel et al 
2001a|) . These methods were developed to overcome the problem of identification of dark matter 
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halos in the very high density environm ents in groups and c lusters. We use a variant of Hierarchical 



Friends-of-Friends Algorithm (HFOF : lKlypin et al.lll999|) for the halo substructure identification. 



To distinguish from HFOF, we call our method HiFOF. The main difference between the HiFOF 
and the HFOF is the way to determine a stable subhalo. The HFOF uses a particle distribution at 
an earlier epoch and checks both the existence and one-to-one correspondence of the progenitor 
particle clusters. A candidate subhalo is considered to be stable if it has one (or two) progenitor(s) 
and if it is the only descendant of the progenitor(s). The HiFOF examines the virial condition of 
subhalos and identifies the virialized subhalos at the highest level as the stable ones. The detail is 
given below. 



The Friends-of-Friends (FOF; see, e.g. JPavis et al.lll985[) algorithm as a base of HiFoF iden 



tifies virialized cluster halos by using a linking length of 0.2 / to link particles as a group if the 

— —1/3 

separation of two particles are shorter than the linking length; here, / = and no is the mean par- 
ticle density in the simulation box. However, FOF is not capable of finding substructures in cluster 
halos. Our algorithm, HiFOF, applies the FOF algorithm with a hierarchical set of linking lengths 
plus a virial condition check on all identified groups on all level. The detail description for HiFOF 
is as follows. We construct sixteen hierarchical levels by decreasing progressively from the linking 
length 0.2 /, the lowest level, to 0.04 1, the highest level. Because of different mass resolution in our 
simulations, we setup up two different criterions for the minimum particle number allowed for a 
group. In the ACDMiooa and the ACDM200, the minimum particle number for a group searched is 
set to be ten, n„,;„ > 10, corresponding 4.1 x 10'" Mq. As in the ACDMiooh and the ACDMiooc, 
n^in is set to be greater than 30, corresponding to 1.9 x 10'" /i"' Mq and 1.55 x 10'° /i"' Mq, re- 
spectively. To examine the resolution effect due to n^m, the analysis in ACDMiooc with n^m = 10, 
corresponding to 5.15 x 10^ /i^' Mq, is additionally included and marked *ACDMiooc- During the 
hierarchical tree construction, any group with particle number less than n^m is discarded. For every 
cluster halo, once a tree of sixteen levels are built, a virial condition check is then performed on 
every member in the tree. The virial condition is calculated by summing up the potential energy 
and kinetic energy in physical coordinates of all particles in a group. A parameter for the virial 
condition check is defined as 



/3 = 2^ - + 1 (1) 



where / sums over the entire particles in a subhalo. It is known that for a virialized system the 
virial parameter yS is equal to 0. However, the virial condition for subhalos within a host halo could 
be made less restrictive due to the potential well of the subhalo residing inside an even deeper host 
halo potential. Nevertheless, it is difficult and impractical to really calculate the extra potential 
energy provided by the host halo. We therefore adopt a loose constraint that yS has to be greater 
than -3; that is, we allow a particle group with Ekmetk < 'AE potential to pass the virial condition. The 
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samples adopting this loose constraint turn out to do better in the comparisons with observations 
than the ones adopting the original condition, /3 = 0. For example, in the ACDMiooa at z = 0.3 
with minimum galaxy mass 4.125 x 10^" h'^ Mq, the identified galaxy (see detail in Section [23] ) 
number density for /3 = -3 is 2.01 x 10"^ Mpc'^ and for = is 0.26 x 10"^ Mpc'^. It can 
be seen that the galaxy density with /3 = is much lower than the one with /3 = -3, and a similar 
situation appears in the galaxy mass functions. 

Additionally, unlike the unbound procedure that iteratively removes the unbound particle s 



with the greatest energy until only bound particles remain in the subhalo (|Springel et al.ll2001a|) . 
our subhalos are treated as a whole and can not be separated into bound and unbound particles. In 
other words, if a subhalo can not pass the virial condition check, we will drop it entirely. After the 
virial conditions of all subhalos in the tree are checked, the virialized subhalos on every branch 
at the highest possible level are then selected to be the samples used in our analysis; other un- 
virialized subhalos and virialized subhalos not at the highest possible level are dropped. Figure [U 
schematically shows how the HiFOF algorithm works. Continuing the same procedure on all FOF 
halos, subhalo samples are then constructed. This full subhalo sample is also called the HiFOF 
sample in this study. 

Figure [2] shows the cumulative mass functions of (bright and dark) subhalos for the four simu- 
lations at redshifts z = (left) at z = 1 (right). The densities are 7. 1 1 x 10"^ Mpc~^ (ACDMiooa), 
7.33 X 10-2 Mpc-^ (ACDM200), 1.53 x 10"^ Mpc-^ (ACDMiooh), 1-53 x 10"' Mpc'^ 
(ACDMiooc), and 4.37 x 10"^ Mpc~^ CACDMiooc) at z = 0. The accumulated subhalo mass 
functions in our simulations basically are consistent with each other. In other words, at the same 
redshift the mass function is mildly dependent on the box size and the mass resolution. However, 
it is observed that the rapid increase appears at the low mass end and is possibly due to the res- 
olution effect when different resolution runs are compared. At z = 1, the four simulations agree 
with each other to a lesser degree. It indicates that the difference likely arises from the different 
softening-length adopted as well as from sample variance. 



2.3. Galaxy Model 

Galaxy formation involves complicated processes. Despite that, cooling is essential to lower 
the specific entropy in the gas, thereby increasing the gas density. Bremsstrahlung and line coolings 



are the most efficient cooling mechanisms for this purpose (e.g., iKaze et al.lll996|) . However, these 
cooling mechanisms require the gas temperature to be above or close to the ionization temperature. 
When a self-gravitating gas is cooled, it can also be heated by adiabatic contraction. It is thus 
possible for the gas to maintain a temperature above the threshold cooling temperature during the 
contraction, thereby yielding to cooling runaway to form stars. In our galaxy model, we hence 
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demand that the host halo potential must be sufficiently deep for the gas to get above the ionization 
temperature. This condition is translated to a threshold host halo mass, only above which stars can 
form within the host halo. 

Now, consider a halo with a virial mass M and a virial radius R, M can be related to R as 

M = —p.,rR' = —Ac(z)pbgR\ (2) 

where pv„- is the virial density. Adz) the mean overdensity of spherically virialized objects formed 
at redshift z, and phg the background densit y = Pn(l + z)^. The over density Adz) at each redshift 
can be evaluated using a fitting formula by Kitayama & Suto ( 19961) . For example, » 335 at 
z = and A, « 201 at z = 1 for our cosmology, Q.o = 0.3 and Q./^ = 0.7. 

In order to facilitate efficient gas cooling, a deep gravitational potential is needed for the virial 
temperature to exceed roughly the ionization temperature, i.e. fcr,,,> ^ ja^m^c^, where a is the fine 
structure constant, the electron mass, and c the speed of light, and b has is a fudge factor of 
order unity. The virial temperature for a star-forming halo can be obtained as 

GMnip b 2 2 

kTyir « — > -a nigC (3) 

5R 2 

where G is the gravitational constant, k the Boltzmann constant, and the proton mass. Using 
Equation (|2]) to replace R by M, it follows that 

\Adz)ptgy"Gm/ 

« 7.0 X 10^2 (A,(z)(l + z)')'"^ h~' Mo, (4) 

where a value of Z? « 0.8 has been adopted. At z = 0, the threshold halo mass Mc « 4.0 x 
10'^ h'^ Mq. We apply the mass threshold to select those FOF halos, which host our galaxy 
samples. All other subhalos hosted by less massive FOF halos are considered dark, void of star 
formation. 



3. Simulation Results and Observations 

In this section we will compare our galaxy samples with the observed galaxies. Comparisons 
include the differential mass function, the galaxy number density evolution, the two point correla- 
tion function at low and high redshifts, the halo occupation distribution, and finally the kinetic pair 
fraction. 
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3.1. Differential Mass Function 

Luminosity function (LF for short) 0(L) gives the relative numbers of galaxies of different 
luminosities, and is so defined that (p{L)dL is the number of galaxies in the luminosity interval 
L ^ L + dL per unit volume of the Universe. The luminosity function of galaxies can be fitted by 
Schechter's (1976) formula, 

= (l>*(^rexpA^. (5) 

J^-k J^-k l^k 

We adopt (P* = {2639t\f^) x lO'^ Mpc-\ a = -1.30, and = iq(~(m*^m^b )I2.5)i^ ^ ^ 937 ^ 



0.25 X lO^^Lo for = -21.07 + 0.13 from DEEP2 B-band galaxies at z ~ 0.3 (iFaber et al.ll2007l) 
for comparison later in our analysis. The mass function and luminosity function are related by 
F{M)dM = (f>{L)dL. If L can be expressed as a function of M when there is no sca tter between M 



and L, the mass function F{M) can be derived from the luminosity function 0(L). iHoekstra et al 



(|2005|) measured the weak-lensing signal as a function of rest-frame B-, V-, and R-band luminosi- 
ties for a sample of "isolated" galaxies from the Red-Sequence Cluster Survey with photometric 
redshifts 0.2 < z < 0.4. They fit the measurements with a power-law for the mass-to-light ratio 

where A//,d is the virial mass of a fiducial galaxy of luminosity L = 10'" h'^L^ Q, and x indicates 
the relevant filter. In B band, they obtained M/,rf = 9.9!l^x 10" h'^ Mq and/3 = 1.5 + 0.3. With 
the power law form of the mass-to-light ratio, F{M) can be found as follows, 

0*(1O'°)"+' M a.i , 

F(M) = "^11 X (—)--' X exp ^-/-^ . (7) 

However, the galaxy mass measured by the weak lensing signal may in fact not represent the 
true galaxy mass correctly. Notably, the weak lensing mass may include the mass of the host halo 
(FOP halo) of an isolated galaxy. We select those from our sample galaxies with only a single 
HiFOF subhalo residing in an FOF halo in our simulations, i.e. the isolated galaxy, to suit the 
observation requirement. These subhalos are a small population of the entire galaxy sample and 
are used as a mass calibrator. The ratios of the selected isolated galaxies to the whole galaxy 
sample are 32% for ACDMuna, 16% for ACDMum, 27% for ACDM2Q0, and 12% for ACDMwoc- 
The result is shown in Figure [3] and the 1 cr error bars are plotted only for ACDMioo^- It is found 
that the relation reveals a power-law form. The power law form follows 

logioiMg) = a* lohioiMi,) + b, (8) 



where a = \ .00621 and b = -0.806326 for the average on all simulations. Extrapolating this power 
law relation to smaller galaxy mass subhalos, we apply this relation to our entire galaxy samples. 
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Figure|4] shows differential mass functions derived from the luminosity function and the M/L 
relation in Equation (|7]) and our galaxy samples at z = 0.3 after applying the Mg - power law 
relation in Figure [3] to our galaxy samples. The mass functions of four simulations agree with 
the observation, and also with each other, quite well. In highest mass range, M > 10^^ Mq, 
and in the low mass end, the profiles of our mass functions show the excess compared with the 
observational data. 

In order to understand the low mass excess, a galaxy mass cut with n„„„ = 10 denoted as 
*ACDMiQQc, in contrast to n„„„ = 30 for ACDMiooc, is analyzed. We find that the low-mass excess 
is pushed toward lower mass. Hence, it is expected that the low-mass excess can be pushed to the 
very low mass end if the simulation mass resolution approaches infinitely high. As for the high- 
mass excess, it is likely related to the over-abundance of CD galaxies at the cluster centers, whose 
population deviates from the Schechter's function. 



3.2. Galaxy Density Evolution 

As the high redshift data are gradually gathered in recent years, study of galaxy properties in 
time evolution becomes feasible. An aspect to test our galaxy model is to make comparison with 
the observed galaxy density evolution. When we integrate Equation (|5]) over luminosity, it gives 
the galaxy number density, 

"(>L) = (^iL')dL'. (9) 
That is, with gi ven 0*, a, L^,, and a l uminosity cut L, the galaxy number density can be obtained. It 



is suggested by iConroy et al.l (120051) that isolated galaxies at z ~ 1 have a similar mass as isolated 
galaxies that are 1 mag fainter at z ~ 0. They found that there has been little or no evolution in the 
halo mass of isolated galaxies with magnitudes in the range ~ M* - 0.5 to - 1.5, even though M* 
has evolved by ~ 1 mag over this redshift range. This result is adopted in our analysis. However, 



Conroy et al.l (120051) also assume that there is no evolution for isolated galaxies with magnitudes 
~ -I- 2.5 logjQ 4, which is equivalent to The luminosity cut L in Equation ^ for n(>^ is 

set to be greater than L^{z)/4-. That is, the cut is L^(z) -dependent. In contrast, we need to find the 
corresponding redshift-independent mass-cut from the simulations to make a correct comparison 
if the previous result and the assumption are to be adopted. Unlike others using the galaxy density 
as an indicator to provide a certain absolute magnitude threshold so as to determine the galaxy 
mass, we take a different approach. We combine the LF of the DEEP2 B-band galaxies and the 
mass-to-light relation at z ~ 0.3, as well as the Mg - A//, relation discussed in Section [XT] to obtain 
the corresponding mass-cut. The luminosity cut of L*(z ~ 0.3)/4 given from the observed DEEP2 



lOOflj 



LF is then converted to a galaxy mass-cut in our samples with 7.02 x\0 h Mq in the ACDM 
7.06 X 101° h-^ Mq in the ACDM200, 7.01 x 10^° h'^ Mq in the ACDMioob , 6.66 x 10"' h'^ Mq in the 
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ACDMiooc, and 5.79 x 10"' h ' Mq in the *ACDMioof • The mass cut is applied to all redshift and 
the redshift evolution of the galaxy number density can then be obtained for our galaxy samples. 

Figure |5] presents the number density evolution of our da ta and the o b serve d galaxy density 



data by integrating the lumi nosity functions ( Equati on |9l) of 



Faberetal 



measurements of Combo- 17 (FaberetaL 2007). PDF (Gabasch et al. 



2003h. VVDS ( Ilbert et al.ll2005h . 2df dNorberg et al.ll2002h . DEEP2 dpaber et al.ll2007h . and SDSS 



(20071), which include 



2004), Bell SDSS (Bell et al. 



( Blanton et al. 20031) . The profiles of our galaxy density are basically in broad agreement with 



the observed galaxies. Our data show the same evolutionary trend, a decline with the redshift. 
However, at low redshift, a depletion is found in our simulations. The depletion can probably be 
attributed to the failure of the assumption for the redshift-independent mass-cut at low redshift, 
and can be corrected if a different mass-cut at low redshift is assumed. 



3.3. Two-Point Correlation Function 

The two-point correlation function ^(r) (CF for short) is the most used indicator to quantify 
the degree of clustering in a galaxy sample. It is defined as a measure of the excess probability 
above Poisson for finding an objects in a volume element dV at a separation r from an otherwise 
randomly chosen object, 

dP = n[l + ^(r)]dV, (10) 
where n is the mean number density of the object in question. It follows a simple power law form, 

^ir) = i-r, (11) 

where ro is the correlation length and y is the power index of CF. 

To learn about the real- space correlation function, we follow the standard practice and com- 
pute the projected correlation function 

Xf max I 
^(^rl + y^)dy. (12) 

The integration limit r,„ax is set to be 40 h'^ Mpc for z = as in SDSS samples. The projected 
correlation functions of our galaxy samples at z = are shown in Figured The so lid line represents 



the volume limited sample {M,. < -19.0) of SDSS galaxy (IZehavi et al.ll2005|) with the density 
1.507 X 10-2 Mpc-\ ro = 4.56 + 0.23 h-^ Mpc, and y = 1.89 + 0.03. We select the mass cuts for 
our galaxy samples to match the number density of SDSS. It is seen that both the amplitude and 
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shape of the projected CF at z = are in very good agreement with those of the SDSS data for all 
simulations. The close agreement of galaxy correlation functions implies that the overall clustering 
of the galaxy population is determined by the distribution of their dark matter subhalos subject to 
the condition of a sufficiently deep halo potential to trigger galaxy formation. In Figure |7] we plot 
the galaxy-galaxy CFs ^ir)gg and the subhalo-subhalo CFs ^(r)/,/, without imposing the condition 
of a sufficiently deep host halo potential for comparison. It clearly shows that ^(r)/,/, also obeys a 
power law form, but with a flatter slope and a smaller amplitude. Note that the number density of 
the subhalo sample that gives ^(r)hh is roughly 5 times higher than that of our galaxy subhalos for 
all simulations. 

Figure [8] shows the results of the projected correlation functions at z = 1 of our galaxy samples 
for the four simu lations and the volume-limited sample of bright, {Mb > -19.0), DEEP2 galaxies 
Jcoil et al.ll2006L . The upper integration limit is set to r,„ax = 20 h ^ Mpc in our simulations to 
agree with the DEEP2 samples. The DEEP2 galaxy sample has a density of 1.3 x 10"^ hr' Mpc~^ . 
We therefore select our samples to match the density of the DEEP2 sample. In Figure [H it can be 
seen that the correlation functions of the DEEP2 galaxies and our samples also agree very well. 

In Figure H] and Figure [lUl the correlation length and the power index y are plotted as a 
function of diff'erent space density at z = and z = 1, respectively. The CFs are fit over the 
range of scales from 0.2 to 13 Mpc where the errors in ^(r) in our sam ples are the "jackknife " 
1 cr errors, computed using the eight octants of the simulation cube ( see [Weinberg et al.ll2004|) . 



Figure [9] is basical ly adopted from Fig ure 11 in Kravtsov et al. ( 2004) which i 
Degre e Field (2dF: lNorberg et allbooj) . SDSS galaxv survevs feehavi etaDl2002 



udes the Two- 



Budavari et al 



20031) and their simulation data. We additionally add the newest SDSS data Jzehavi et allboosl) 



and our galaxy data into Figure |9l The strong dependence of correlation length ro on the number 
density n is evident in Figure[9l That is, brighter galaxies are more clustered. In addition, a nearly 
constant profile of y in observations is seen. The amplitude of y in our simulations, in general, 
also presents a roughly flat trend although din ^ 10"^ Mpc'^ our n values sho w slight depletio n 
compared with the obse rvations. At z = 1 we plot the DEEP2 data released by lCoil et al.l (|2006h . 
the simulation results of iKravtsov et al.l (|2004h . and our data. A similar conclusion as at z = is 
obtained for ro, and our values of y also show good consistency with the DEEP2 galaxies. 



3.4. Halo Occupation Distribution (HOD) 



The HOD formalism, developed during the last several years, has become a powerful theo- 
retical framework for predicting and interpreting galaxy clustering. The original HOD uses the 
probability P(N\M) to describe the b i as of a class of galaxies that a halo of virial mass M con- 



tains N such galaxies. iBerlind et al.l (|2003|) studied HOD and found that at a given halo mass it 
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is statistically independent of the halos large-scale environment. In addition, they compared HOD 
of a semi-analytic model and of gasdynamics simulations in detail and concluded that the semi- 
analytical HOD for samples of the same space density agree remarkably well with simulations, 
despite that the two methods predict different galaxy mass functions. 



Kravtsov et al.l (|2004) showed that HOD can actually be understood as a combination of 
the probability for a halo of mass M to host a central galaxy and the probability to host a num- 
ber of satellite galaxies that obeys Poisson statistics. We analyze the first moment of HOD, 
{N(M,,)) = (NciM,,)) + (A^,(M/,)), as a function of host mass M/, for the halo samples, where (Nc) 
is the HOD of central galaxies modeled as a step function {A^^.) = 1 for M > M„„„ and (Nc) = 
otherwise , and (A^,) is number of " the satellite galaxies modeled as a power law (A^.,) oc M" sug- 
gested by iKravtsov et al.l (|2004h . IZehavi et al.l (|2005|) used this HOD framework to interpret their 
results. Due to a finite (but important) subset of information encoded in the correlation function for 
galaxy clustering, a restricted HOD model was employed and fitted with a small number of free pa- 
rameters. The HOD formulation they implemented has three free parameters: M^m, the minimum 
halo mass for galaxies above the luminosity threshold. Mi, the mass of a halo hosting one satellite 
galaxy above the luminosity threshold, and a, the power law slope of the satellite mean occupation 
function. They took M„„„ to be fixed by matching the observed space density of the sampl e, leaving 
Ml and a as free parameters to fit their projected CFs. More recently Jzheng et aL ( 2007 ) modeled 
the luminosity-dependent projected two-point correlation function of DEEP2 (z ~ 1) and SDSS 
(z ~ 0) galaxies within the HOD framework. They adopted a more flexibl e param eterizations with 
five parameters, motivated by the less satisfactory results of IZheng et al.l (|2005|) as well as of the 
three-parameter parameterizations. In the following test on this aspect, only the three-parameter 
parametrization scheme is adopted for simplicity. 

We compute the HOD of our galaxy samples to obtain M„„„, Mi, and a as a function of galaxy 
number density to make comparisons with SDSS and DEEP2 galaxies. In Figure [TTT a'). from top 
to bottom, our results of M„„„, Mi, and a in the four simulations at z = and those of SDSS 
galaxies are compared. It can be seen that M^m and Mi of four simulations not only agree with 
each other, implying the lesser degree of dependence on the box size and the resolution, but also 
match the observation quite well no matter in the magnitude or in the trend. Despite the increasing 
trend of the slope a with a decreasing density is similar to that of SDSS galaxies, the a value of our 
HODs on the whole appears to be slightly smaller than that of the observations. This may results 
in slightly insufficient subhalos found in massive clusters in our samples. Figure \TT\ b) shows our 
results at z = 1 and the observational data of DEEP2 galaxies. As on can see, the agreements on 
M^in, Ml, and a are also good. 
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3.5. Pair Fraction 



Since our galaxy sample can well reproduce the observed evolution of galaxy number density, 
the mass function, the 2-point correlation function and the HOD, we therefore make an attempt to 
put another comparison to the simulation resolution limit by exploring the pair fraction and merger 
rate. 



DEEP2 Team (iLin et al.1 120041) explored the kinematic close pair fraction and the merger 
rate up to redshift z ~ 1.2 and disclosed weak evolution in the galaxy pair fraction. Assuming 
mild luminosity evolution, they found the number of companions per lumi nous galaxy to e volve 
as (1 + zf, with m = 0.51 + 0.28 for the r„^„ = 50 kpc case. Recently, iLin et all (l2008b used 
more complete data to study the same problem, and obtained an improved result with m = 0.41 + 
0.20 for all galaxies. These two studies consistently reveal that the pair fraction of galaxies indeed 
undergo weak evolution. 

We study the pair fraction problem with our galax y samples; our appr oach is in variance with 
applying a "hybrid" formalism to address this problem (iBerrier et al. which combines an N- 

body simulation to account f or large-scale structu re and the host dark matter halo population and an 
analytic substructure model (|Zentner et al. to identify sa tellit e galaxies within the host halos. 
We test our samples with the same setup as in Lin et al. (2004) and Lin et al. (2008). They defined 
close pairs such that the projected separations satisfy 10 kpc < Ar < r^^x, where r,„ax = 30, 50, 
or 100 h'^ kpc, and the rest-frame relative velocity Av less than 500 km s^^ . To ensure the selection 
of the same types of galaxies at different redshifts in the presence of luminosity evolution, a specific 
range in the evolution-corrected absolute magnitude M^, defined as Mb + Qz, was adopt ed in our 
analysis where the evolution is parame terized as M{z) = M(z = 0) - Qz. iLin et al.l (|2008h adopted 
2 = 1.3 found by iFaber et al.l (|2007|) and restricted their analysis to galaxies with luminosities 



-21 < < -19 for z = 0.45-1.2. Applying the mass-to-light ratio in Equation ^ for B-band 



with Mfid 



9.9t\ix 10^1 ri Mg 



and /3 = 1.5 + 0.3, we are able to convert the mass to the 



luminosity. The luminosi ty range, -21.39 < M n < -19.39, at z = 0.3 is considered because the 
mass-to-light ratio data of iHoekstra et al.l (|2005[) was collected at an average redshift z ~ 0.3. The 
mass range at z = 0.3 after conversion is found to be 2.38 x 10" h'^ Mq< M < 3.78 x 10'^ Mq. 
Moreover, the average number of companions per galaxy is defined as 

2N„ 



(13) 



where Np is the number density of ind ividual paire d galaxies and the number density of galaxies 
as suggested in iBerrier et al.l (|2006f ). Following iBerrier et al.l (|2006l) . we also adopt a fiducial 
relative line-of-sight velocity difference of |Av| = 500 km s~^. A close-pair cylinder volume is then 
constructed by a constant Av cut. In our analysis, only the r^ax = 50, and 100 kpc cases are 
considered due to the limitation of our mass resolution. 
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Figure [12] shows the redshift z versus the pair fraction for Vmar = 50 (bottom) and r, 



m 



2003 



2003 



Driver et al. 



100 /i ^ kpc (top ). The pair fracti on data are taken from lLin et al.1 (l2004h incl uding SSRS2 (Ida Costa et al. 
19981). CN0C2 dYee et alJ boOOl). and DEEP2 early data doavis et al.ll2003l) and from full samples 



Lin et all (12008 ) including SRSS2, CN0 C2, MG C (Millennium Ga laxy Catalog, iLiske et al 



20051 : kUen et al.ll2006l) . TKRS (Iwirth et al.ll20o3 ). and DEEP2 JPavis et al. 



20071). Our data points are consistent with the DEEP2 medium to high redshift observations 



although at low redshift we have lower values than other observatio ns. Howeyer, the mild evolution 



trend i s seen. Using the fitting formula, A^c(0)(l + z)'", proposed in iLin et al.l (|2004t) and lLin et al. 
(|2008|) . we find the following results. For Vmax = 50 h'^ kpc , NdO) = 0.036 and m = 1.17 in 
the ACDMma , NM = 0.074 and m = 1.07 in the ACDM20Q, NM = 0.046 and m = 0.96 in 
the ACDMum, NdO) = 0.025 and m = 1.43 in the ACDMiooc, and NdO) = 0.028 and m = 1.37 
in the *ACDMwoc- For r,„ax = 100 h-^ kpc, NdO) = 0.137 and m = 0.65 in the ACDM^oa, 
NdO) = 0.207 and m = 0.95 in the ACDM200, NdO) = 0.182 and m = 0.48 in the ACDMimt, 
NdO) = 0.107 andm = 0.98 in the ACDMiooo and A^c(O) = 0.1 16 a nd m = 0.92 in th e *A CDM^mr 



Appar ently, our pair fractions show stronger evolution than those in lLin et al.l (|2004|) and lLin et al. 



(120081), where m = 0.51 + 0.28 and m = 0.41 + 0.14 for v,,^ = 50 h'^ kpc, and m = 0.47 + 0.18 



and m = 0.29 + 0.05 for r^ax = 100 kpc, but appear to agree with the results of iBerrier et al. 
(120061) for Vmax = 50 h'^ kpc, where m = 0.42 + 0.17 for Vnow, the maximum circular velocity of 
the subhalo at the current epoch, and m = 0.99 + 0. 14 for Vm, the maximum circular velocity of the 
subhalo when it was first accreted into the host halo. 

The steeper slope in our sample is rooted on the too-low galaxy density at low z, which can 
basically be attributed to the failure of our assumption on the redshift-independent mass-cut. We 
note that the pair fraction critically depends on the galaxy density and the density is determined by 
the luminosity cut. Consequently, if a wrong luminosity evolution assumption is adopted, this can 
lead to a high galaxy mass-cut, a lower galaxy density, and finally a lower pair fraction. Our data 
at low redshift reveal this effect and result in a steeper slope. 



4. Conclusion and Discussion 

In this study we have used the HiFOF method to locate and identify subhalos, and we advance 
the idea that galaxies form only in a sufficiently massive halo due to its sufficiently high gas cooling 
efficiency. It is estimated that the threshold halo mass Mc ~ 4.0 x 10'^ H'^Mq at z = 0, and 
our galaxy samples are constructed from the HiFOF subhalos embedded in the halos satisfying 
this mass threshold. There are far more subhalos that are hosted by halos not meeting the mass 
threshold. They are dark halos where star formation cannot occur. To test our model, we have 
analyzed the differential mass function, two-point correlation function (CF), the space density 
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evolution, the HOD, and the pair fraction from our galaxy samples. The results are summarized as 
follows. 

(1) Combing the luminosity function (LF) and the mass to light relation (M/L), a differential 
mass function can be derived. Based on the lensing signal from Red-Sequence Cluster Survey to 
determine the galaxy mass, we have corrected the lensing mass in our galaxy samples so as to 
make fair comparison. We find that the mass functions of our four simulations agree with the LF 
and M/L data very well. 

(2) We compare our galaxy density evolution with the evolution of the observational galaxy 
number density through the integration of a LF with a proper luminosity cut. Our simulations 
basically show fair consistency with observations, except at low-redshift. The insufficient density 
is likely due to the failure of the assumption at low redshift about the redshift-independent mass- 
cut. The depletion in the galaxy number density also reflects a similar result in pair fractions found 
at low-redshift. 

(3) We have analyzed CFs of our galaxy subhalos at z = and at z = 1. At z = 0, we 
select our galaxies to have the density 1.507 x 10"^ Mpc'^ to make comparison with the volume 
limited sample {Mr < -19) of SDSS galaxies. Both the slope and the amplitude are consistent 
with the SDSS. At z = 1, we select samples at n = 1.3 x 10"^ Mpc~^ to setup the same density 
condition as in the volume-limited sample of bright (corrected), Mb > -19.0, DEEP2 galaxies. 
The projected correlation functions are compared, and our profiles are similar to the one found in 
DEEP2 galaxies. 

(4) We also study ro and y as a function of the space density n. Our galaxies reveal the 
increase of ro as n decreases, a similar result also found in the observation data. In other words, 
the correlation length ro has luminosity dependence, and the brighter the galaxies, the larger the 
correlation length. On the other hand, the power index y of observation data displays a nearly 
constant behavior over a wide range of n at z = 0. Despite a smaller values of y are obtained 
around n x 10'^ in our simulations, a roughly constant trend on the whole is seen. The same 
results are reached at z = 1 as well. 

(5) We have parameterized our HODs with three parameters Mi, M,„in, and a as a function 
of the galaxy number density. Our values of Mi and A/,„,>, and their increasing trend as the galaxy 



density decreases show good agree ment with SDSS data (|Zehavi et al.ll2005l : IZheng et al.ll2007|) 



and DEEP2 data (IZheng et al.ll2007|) . and the depletion of the slope or at z = may be attributed to 



insufficient subhalos found in massive clusters. 

(6) We have analyzed the kinematic pair fraction of our galaxy samples by selecting the mass 
in between 2.28 x 10' ' Mq and 3.61 x 10'^ Mq which are obtained by applying the observed 
mass-to-light ratio to convert the luminosity range -21.39 < Mg < -19.39 at z = 0.3. Parame- 
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terizing the evolution of the pair fraction as (1 + z)™, iLin et alj (|2008|) find that when r^ax = 50 
kpc, m = 0.41 + 0.14 for the full sample. We obtain m = 1.17 in ACDMiooa, m = 1.07 in 
ACDM200, m = 0.96 in ACDMioot, m= 1.43 in ACDM^ mr, and m = 1.37 in *ACDMiooc- When 
fmax = 100 kpc, in contrast to 0.29 + 0.05 found by iLin et al.l (|2008|) . we find m = 0.65 in 
ACDMwQa, m = 0.95 in ACDM200, m = 0.48 in ACDMimb, m = 0.98 in ACDMiooc and m = 0.92 
in *ACDMiooc- The steeper slope basically results from insufficient close pairs found at low red- 
shift, which is interpreted as the deficient low-z galaxy density in our samples due to the failure of 
the assumption on the redshift-independent mass-cut. However, the situation is less severe for the 
fmax = 100 kpc case where milder evolution trend is more evident. 

Most of the above tests for our galaxy samples provide supporting evidence for our galaxy 
model. However, this model has some problems of its own. First, this model do es not provide 



inform ation of the galaxy morphology. Second, unlike the Millennium simulation (jSpringel et al 



20051), this model does not contain any empirical ingredient to trace individual galaxy evolution. 



The solution to the first problem requires the gas component to be included in the simu- 
lation, which can reveal the amount of residual gas angular momentum in the subhalo thereby 
providing additional information of morphological types. The complexity of the second problem 
is intractable from the first-principle calculations and it needs an empirical model, much as what 
was added to the Millennium, to approximate the physics involved. 

We would like to thank Lihwai Lin and Tak-Pong Woo for many helpful discussions and 
suggestions. This project is supported in part by the grant: NSC-92-2628-M-002-008-MY3 (1/3). 
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Fig. 1 . — Example of the HiFOF algorithm is illustrated with a five-level hierarchical tree structure 
of one FOF halo. Points with "V" stand for "virialized" halos whereas points with "N" are for "un- 
virialized" halos. Blue, brown, and light-blue points represent the FOF halo, the selected galactic 
subhalos, and the subhalos dropped. As mentioned in the text, the hierarchical tree is constructed 
through the FOF method with gradually reduced linking lengths. To determine whether a subhalo 
is bound, the virial condition of it is checked. We start with selecting the samples from virialized 
subhalos on the highest level; once a member in a branch is selected, others in that branch will 
be discarded. We then move to the next lower level and repeat the procedures just described. A 
completed subhalo sample is then acquired after every level is examined. 
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Fig. 2. — Cumulative mass functions of (bright and dark) subhalo samples obtained in our anal- 
ysis at redshift z = (left panel) in the ACDMiooa (solid line), ACDM200 (long-dashed line) , 
ACDMioob (short-dashed line) , and ACDMiooc (dotted line) and at z = 1 (right panel). The mass 
range plotted are down to the minimum mass of 4.1 x 10'" /z ' Mq in the ACDMu^a and ACDM200, 
to that of 1.9 X 10'" h-^ Mq in the hCDMm)h, to that of 1.52 x 10'" Mq in the ACDMunc, and 
to that of 5.16 X 10^ Mq in the *ACZ)Miooc- These mass functions are generally consistent with 
each other. As a result, we believe that the influence of the box size and resolution effect are not 
important for our HiFOF samples. 
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10^^ 10^^ 10' 



Fig. 3. — Diagram shows Mf, (Mass of the isolated FOF halos) vs. Mg (Mass of the HiFOF 
identified subhalos) at redshift z ~ 0.3. As shown in the figure, the relation can be approximated 
as a power law. This power law is used to provide the correction on galaxy mass from the weak 
lensing mass. The error bars represent 1 cr values, and similar bars also exist for other curves, but 
omitted for the clarity of presentation. 
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Fig. 4. — Differential mass functions at redshift z = 0.3 of LF + M/L data (solid line), the 
ACDMiQQa samples (long-dashed line), ACDMjqo (short-dashed line) , ACDMioot (dotted line), 
ACDMiooc dashed-dotted line) , and *ACDMioof (short-dashed dotted line). The mass of our 
galaxy samples is modified according to the relation in Fig|3l We list the space density as follows. 
The number density n = 2.01 x 10"^ Mpc'^ in ACDM^oa, 2.1 1 x 10"^ Mpc'^ in ACDM200, 
3.22 X 10-2 Mpc-^ in ACDMym,, 2.82 x 10"^ Mpc-^ in ACDMmc, and 4.77 x lO'^ Mpc'^ 
in * ACDMiQQc- Our data profiles are similar to the observational data line despite deviations shown 
in the two ends, the lower and higher mass range. 
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Fig. 5. — Galaxy den sity as a function o f the redshift is plotted. Observational data are based on 
the data of the LFs in iFaber et al.l (|2007[) and can be converted to the galaxy density evolution by 
using Equation A proper luminosity cut, L > L*(z)/4, is adopted for the density calculation, 
and combining the data of the DEEP2 B-band galaxies and the mass-to-light relation at z ~ 0.3, as 
well as the Mg - A//, relation, a corresponding mass-cut can be obtained and then applied to galaxy 
samples at other redshifts to obtain n{z). 
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Fig. 6. — Projected correlation functions of the galaxy samples at z = in the ACDMiooa 
(dashed line), ACDM2Q0 (dotted line), ACDMioo/, (open triangle with the error bar), and ACDMiooc 
(dashed-dotted line) are compared with the proj ected correlation fu nction of the volume lim- 
ited sample (Mr < -19) of the SDSS galaxy Jzehavi et all llOOSh with the density 1.507 x 
10-2 Mpc'\ 




Fig. 7. — The correlation functions ^{r)gg of the galaxy samples at z = in the ACDMi^oa (solid 
line), ACDM200 (dashed line), ACDMiooi (short-dashed line), and ACDMiooc (dotted line) are 
plotted. In contrast, the subhalo-subhalo projected correlation ^(r)/,/, is presented with #. It is 
clearly shown that ^(r)hh obeys a power law form, but with a flatter slope and a smaller amplitude. 
That is, our galaxy samples behaves substantially differently from HiFOF samples. 
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Fig. 8. — Comparison of the projected correlation functions between our four simulat ions at z = 1 
and the DEEP2 volume limited data of bright, Mb < -19.0, galaxies JCoil et allboOQ) . The CF of 
the DEEP2 data has a density 0.013 Mpc~^ . We therefore select our samples to have the same 
density, and the results can be found to be in good agreement with the DEEP2 data. 
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Fig. 9. — Top: Best-fit correlation scale ro as a function of number density at the present-day 
epoch. Bottom: Best-fit slope y of the correlation func tion as a function of number density. This 



plot i s based on the Figure 1 1 in iKravtsov et al.l (|2004|) plus the newest SDSS data dZehavi et al 



12005 ) and our data. Our values of the correlation length ro are consistent with observational data. 
The values of the power index y are lower at density around n = 10"^ Mpc~^. However, the 
agreement between observations and our data is seen. 
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Fig. 10.— Sam e as Figure |9l but for z = 1 . We additionally include DEEP2 data released by 



(|Coil et al.ll2006h as well as our data. It is found that the correlation length and the power index in 
our simulations match DEEP2 data well. 
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Fig. 11. — From top to bottom, the three fitting parameters Mi, M,„jn, and a of the HOD formula 
tion are plot ted as a function of the galaxy number density ng. In [TIT a), the open circ les represent 



SPSS data (|Zehavi et al.l 120051) and the stars with errorbars give also SPSS d ata (|Zheng et al 



20071), and in fTTT b) the stars with errorbars give PEEP2 data dZheng et al.ll2007h . while the filled 
square, circle, triangle, diamond, and pentagon are for the ACDMiooa, ACDM200, ACDMioofc, 
ACDMiooc , and *ACDMiooc data. Consistency between SPSS galaxies and our galaxy samples is 
observed at z = 0. The good agreements are also seen between PEEP2 galaxies and ours at z = 1. 
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Fig. 12. — Pair fraction, Nc, as a function of redshift z. The crosses mark the results from DEEP2 
early data (2003). The open squares represents measurements from SSRS2, MGC, and CN0C2. 
The open circles show DEEP2 data (2007). The open triangles are from TKRS. The stars, filled 
squares, circles, triangles, and diamonds give our simulation results from ACDMiooa, ACDM200, 
ACDMum, ACDMiooc, and *ACDMiooc. Two cases of r^^, 50 (bottom) and 100 (top) kpc, 
are considered. 
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Table 1: Simulation Parameters 





Lbox 






fz=0 




Name 


(h-^ Mpc) 




ih-' Mo) 


( h-^ kpc) 




ACDMma 


100 


256^ 


4.125 xlO^ 


6.0 


0.94 


ACDM200 


200 


512^ 


4.125 xlO^ 


10.0 


0.94 


ACDM.oob 


100 


512^ 


6.188x10^ 


3.0 


0.94 


ACDMiooc 


100 


512^ 


5.156x10^ 


10.0 


0.94 



